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We describe a systematic expansion for full QCD. The leading term in the expansion gives the valence approx- 
imation. The expansion reproduces full QCD if an infinite number of higher terms are included. 
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1. INTRODUCTION 

The observation that recent valence approxi- 
mation predictions are quite close to experiment 
strongly suggests that, at least for some observ- 
ables, the systematic error arising from this ap- 
proximation is quite small. It then seems natural 
to ask vifhether an independent quantitative es- 
timate can be made of the systematic error of 
the valence approximation. In this article, we de- 
scribe an expansion for full QCD ||l| which has the 
valence approximation as its leading term, repro- 
duces full QCD exactly if an infinite number of 
higher terms are included, allows the numerical 
calculation of lattice QCD quantities for any fi- 
nite truncation of the expansion, and allows the 
estimation of the error induced by any such finite 
truncation. 

The basic idea behind the scheme we suggest 
is to develop an effective action to approximate 
the log of the determinant of the fermion hop- 
ping matrix. The valence approximation replaces 
this term by a quantity proportional to the pure 
gauge action. Improved approximations can be 
generated by replacing the log of the determinant 
by a more complicated function of the links. A 
useful set of such functions is generated by sums 
of traces of products of links about closed paths 
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on the lattice, and our approximation scheme in- 
volves expanding the determinant of the hopping 
matrix as a linear combination of these functions. 

2. DEFINITIONS 

We will work with Wilson's formulation of lat- 
tice QCD. If M{k) is Wilson's hopping matrix 
for a single flavor fermion with hopping param- 
eter K, then our aim is to find an approxima- 
tion for BetM (= VDetA/tAf). Numerical sim- 
ulations of QCD work on finite lattices, with k 
not too big. For almost all configurations on 
such lattices M'I'Af is positive definite and has fi- 
nite, strictly positive minimum Amin and maxi- 
mum Ajnax eigenvalues. Also, it is simple to see 
(for example by considering a hopping parameter 
expansion) that Det (iW'^A'/) can be expressed as 
a finite linear combination of sums of traces of 
products of links about closed paths. These two 
observations allow us to define log (Det (^M^M) ) 
as an infinite series, convergent with respect to a 
suitably defined norm [|l| 

oo 

log Det (M^Af) = Tr log M^M = ^ a,5, ([/) (1) 

i=Q 

where the Si{U) are a maximal linearly indepen- 
dent set of functions which involve sums of traces 
of products of hnks about closed paths. 

Crucial to obtaining a useful truncation of the 
expansion in Eq. (^ is the ordering assigned to 
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the sequence of Si{U). We choose to partially 
order these functions by the length of the closed 
paths involved. For what follows we need explicit 
forms only for the constant (path length 0) term 



So{U) = 1 

and the path length 4 term 
Si{U) = i^RcTrf/n 



(2) 



(3) 



involving a sum over all plaquettes of the real 
part of the trace of the product Un of links about 
a plaquette. 

3. APPROXIMATION SCHEME 

The series expansion for Tr log M^M defined in 
(|^) allows us to develop a series of approximations 
for full QCD expectation values of observables F. 
If we include two equal mass flavors of fermions, 
(the case most simply treated with exact algo- 
rithms), we have 

{F) = Z-'^ J d^i F exp{(3Si+TTlogM'iM) 

= Z-^ J dn F exp(^(3Si+J2^^^i) (4) 

We approximate these expectation values by 
truncating the series for TrlogAffM. 

If Ln = X]"=o '^i^i such a truncated series, 
and if i?„ = Trlog Af'I'Af — L„ is the remainder 
after truncation, then we define our approximate 
expectation after this truncation by 



Z-^ J dfiF exp{pSi+Ln) 



(5) 



The systematic error induced by the truncation 
is (F) — (F)„, and if the correlations of F and i?„ 
are small, then an estimate of this error is 

(^^> - (^^>„ « {{F (F)„) (i?„ - J>„ (6) 

Still to be determined at this stage are val- 
ues for the coefficients a^. The norm implicit in 
Eq. d) yields 



hm / Tr log M^Af - ^ a^S'J \ = 0. 



(7) 



We can find approximate values for the coeffi- 
cients at any given truncation by minimizing 



(i?2)„ = ((TrlogMtAf-L„)2)„ 



(8) 



as a function of ap . . . a„. In the limit n oo 
the coefficients will approach their exact values 
determined by (|l|) or equivalently (^. 

4. CALCULATING TrlogAftM 

In order to evaluate the coefficients (by min- 
imizing (H)), or to estimate the systematic er- 
ror by equation (^), we must be able to evaluate 
the expectation value of Tr log MtM over a set of 
link configurations. We use two different ideas to 
evaluate these expectations. Firstly, we estimate 
TrlogMtAf as 



Tr log M^M =—J2^l\og Mhl ■ 0, 



(9) 



were (j)a are iV^ pseudofermion vectors each with 
independent random Gaussian components. Sec- 
ondly, to calculate log Af^Af • we use a Cheby- 
shev polynomial iterative scheme 



N 



log M.M . Y-O-nik-"'""' 

+ 5\ogM^M 



/A n 



(10) 



where the T* are Chebyshev polynomials, the hi 
are constants which can be calculated using stan- 
dard methods for orthogonal polynomials (2| , and 
5 is bounded by 



(5< 2exp (-2Ny/X~J\ 



(11) 



Convergence of this series is determined by (5, 
and since log M'^M -(p must be calculated a num- 
ber of times in order to evaluate the trace we 
have found that the most efficient implementa- 
tion of these ideas is first to calculate Amin and 
Amax for each configuration being analyzed us- 
ing a Lanczos algorithm then to evaluate 
dog AftAf •(/) with a Chebyshev polynomial ap- 
proximation whose order N is defined to keep 
-^■\/ Amin /Amax On each configuration constant. 
We have also found that it helps to precondition 
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the log M'I'M calculation by premultiplying M by 
the inverse of a free hopping matrix Afg (kq) which 
has hopping parameter kq, and has all links set 
to the identity. This preconditioner has no ef- 
fect on the physics since, if = M^^{kq)M , we 
have DetA^W = const x Y)etM^M. However, sim- 
ple tuning of kq reduced the work to calculate the 
determinant in our example calculation by almost 
a factor of two. Details of our exact implementa- 
tion of these ideas can be found in . 

5. EXAMPLE APPLICATION 

As a first demonstration of these ideas, we have 
studied a 6"* QCD simulation at /3 = 5.7 with two 
equal mass Wilson fermions of hopping parameter 
K = 0.16. We considered only the first two terms 
in our expansion 

TrlogiVW^Li = ao + ai5i([/). (12) 

Observable expectations in this approximation 
have the form 

{F)^ ^ Z-^ j d^iF exp{{P+ai)Si{U)+ao){13) 

This expectation value is a valence approxima- 
tion calculation but we now obtain an explicit 
number for the shift in (3 compensating for omis- 
sion of the fermion determinant. The coupling 
at which we simulate becomes /3sim — P + ^P, 
where SP = ai is the effective shift in the coupling 
due to the fermion determinant. To calculate Sp 
we generated 160 independent configurations us- 
ing the Cabbibo-Marinari-Okawa algorithm with 
Psim = 5.7. On these configurations, we calcu- 
lated TrlogTVW using 140 independent 0's per 
configuration. 

The major result of all our analysis is a pre- 
diction for the constant 6P which is generated by 
minimizing (|). We found Sp = -0.261 ± 0.014. 
We interpret this result to mean that a valence 
calculation at P — 5.7 is approximately equiv- 
alent to a full QCD calculation at /? = 5.439 
(5.700—0.261) for two equal mass flavors of Wilson 
fermions with hopping parameter k = 0.16. The 
agreement is approximate because we have trun- 
cated the expansion of TrlogiVW at two terms. 

To confirm the approximate equivalence, we 
have also calculated the expectation values of the 



plaquette operator P ~ 1 — iReTr?7n in a va- 
riety of different ways. A pure gauge calcula- 
tion at p = 5.439 gives (P)^ = 0.5190 ± 0.0004, 
A full QCD hybrid Monte Carlo calculation at 
P = 5.439 with two flavors of Wilson fermions 
gives (P) = 0.4451 ± 0.0004. Our approxima- 
tion using the P = 5.7 valence calculation gives 
(P)^ = 0.4499 ± 0.0004 with a systematic error 
of 0.0000 ± 0.0025. The fact that our estimate of 
the systematic error has a central value of 0.0 is 
simply a consequence of our algorithm for find- 
ing sp since we have chosen SP explicitly to make 
this central value equal to 0.0. Other observables 
will have non-zero systematic shifts which can be 
calculated in exactly the same manner. 

6. COMMENT 

The major question which we have not an- 
swered yet is whether our scheme is more efficient 
computationally than simulating with an exact 
full QCD algorithm. In our simple test the extra 
work required to evaluate Tr log M^M turned out 
to be about equivalent to the time required to 
generate a single decorrelated configuration with 
the Cabbibo-Marinari algorithm. How this might 
change as the size of the lattice increases, and 
as the gauge coupling and quark hopping param- 
eters become more physical is still to be deter- 
mined. 
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